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Bosonic two-dimensional self-bound clusters consisting of AT atoms interacting through additive 
van der Waals potentials become unbound at a critical mass mi N ^ ; mi N ^ has been predicted to be 
independent of the size of the system. Furthermore, it has been predicted that the ground state 
energy En of the JV-atom system varies exponentially as the atomic mass approaches m*. This 
paper reports accurate numerical many-body calculations that allow these predictions to be tested. 
We confirm the existence of a universal critical mass m» , and show that the near-threshold behavior 
can only be described properly if a previously neglected term is included. We comment on the 
universality of the energy ratio En+i/Em near threshold. 



PACS numbers: 

Restricting the motion of particles to one or two di- 
mensions can lead to properties that differ dramatically 
from those in three dimensions. The most prominent 
two-dimensional (2D) system is the surface of bulk mat- 
ter. Another example are one- or two-atom layer thin 
films, e.g., atomic or molecular hydrogen films 0, 0], 
grown on substrates. Neglecting the adatom-substrate 
interaction, many properties of such systems can be un- 
derstood within a truly 2D framework. In addition to 
homogeneous 2D systems.it is interesting to consider 2D 
clusters (see, e.g., Refs. |3j |4j). What happens when a 
finite number of atoms is restricted to 2D space? In- 
homogeneous 2D systems can potentially be studied by 
placing a few atoms on the surface of a substrate or by 
confining atoms by external potentials. Effectively 2D 
atom traps have been realized recently |5j, |6J ; extension 
to optical lattices with only a few atoms per lattice site 
is possible with today's technology. These systems are 
particularly interesting since Feshbach resonances allow 
the interaction strengths to be tuned through application 
of magnetic fields. 

Bosonic 2D systems interacting through short-range 
two-body potentials that support one zero angular mo- 
mentum bound state are predicted to exhibit intrigu- 
inguniversal, that is, model- independent, behaviors |7j, 
EM ElElElElHElElElEi i) 2D clus- 
ters with N particles [9J, [UJ, [l2( are predicted to be- 
come unbound when the mass m reaches a critical value 
' A ' • This critical mass is predicted to be universal, 

(AT) 



m 



i.e., 771* = to* [lOj, and to be the same for the cor- 
responding homogeneous system d Il9| . ii) For a given 
system size, the ground state energies En near threshold 
are predicted to change exponentially as the atomic mass 
m decreases 0, 0] . Similarly, for a given atomic mass, 
the ground state energies near threshold are predicted to 
change exponentially with varying system size E3- iii) 
The ratio between the ground state energies of a 2D sys- 
tem with N+l atoms and those of a system with N atoms 
reaches, in the limit of zero-range interactions, a con- 
stant. This constant has been determined analytically for 



small 2D systems: E$/E$ = 16.52 |1E1 and Ef/Ei = 
11.94 18j. For large systems, the ratio E S N+1 /E S N has 
been predicted to approach 8.57 E3- These predictions 
should apply to systems interacting through short-range 
potentials if \En\ <C h 2 /{mr 2 ), where r e denotes the 
maximum of the characteristic length of the two-body 
potential [in our case, the van der Waals length r v( jw, 
r v dw = (mCe/h 2 ) 1 / 4 , where C§ denotes the leading van 
der Waals coefficient] and the absolute value of the effec- 
tive range r c g. 

This paper reports the energetics of self-bound inho- 
mogeneous 2D systems near threshold, and tests under 
which conditions predictions i) to iii) apply. To this end 
we perform numerical many-body calculations, which su- 
persede earlier variational calculations 0, |nj and cover 
a wider parameter range. For N — 2 and 3 we deter- 
mine the ground state energy using basis set expansion- 
type calculations. Since these techniques become un- 
tractable for larger systems, we resort to Monte Carlo 
techniques for N > 3. We show that a proper descrip- 
tion of the ground state energies for small clusters with 
N = 2 — 5 atoms near threshold requires, in addition to 
a term proportional to (to — ml ), a term proportional 
to (to — m^) 2 . For larger droplets with N = 6 and 
7 atoms we only include the linear term; our diffusion 
Monte Carlo (DMC) energies do not allow the term pro- 
portional to (to — to i NS> ) 2 to be determined. We speculate 
that DMC energies covering a larger mass range (see be- 
low) would show that near threshold the quadratic term 
is non-negligible for systems with N > 6. Our study con- 
firms the existence of a universal critical mass to* |l0| . 
This critical mass depends on the two-body potential but 
appears to be universal for bosonic TV-particle systems in- 
teracting additively through a given two-body potential. 
Finally, we investigate the ratio between the total ground 
state energy of 2D systems with N + 1 and N particles. 
Our results are consistent with predictions based on zero- 
range treatments H 0, 0, 0] but do not conclusively 
confirm them for N > 3. We comment on the relevance 
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of the energy scale associated with the effective range. 

Consider the Hamiltonian H for N particles with mass 
m. 



h 2 N A 

2m ^ 1 

i—l i<j 



(i) 



where V 2 denotes the 2D Laplace operator of the zth 
particle, the distance between atom i and j, and V 
the atom-atom potential. We use the realistic helium- 
helium potential without retardation developed by Ko- 
rona et al. [20j (in the following referred to as KORONA 
potential) as well as a Lenard-Jones (LJ) potential with 
parameters e and a chosen to approximate the helium 
dimer interaction, i.e., e = 10.22 K and a = 2.56 A [2l| . 
Both these two-body potentials support a single two- 
body bound state with vanishing angular momentum. In 
the following, we calculate the lowest energy with van- 
ishing angular momentum of the many-body system. We 
first treat "true" bosonic helium systems in 2D, i.e., we 
choose to = mne (here and in the following, He denotes 
the isotope 4 He). To explore the threshold regime, we 
then successively reduce the atomic mass while leaving 
the interaction potential unchanged. 

To start with, we solve the scaled radial Schrodinger 
equation for two atoms by diagonalizing the Hamiltonian 
using B-splines. For each mass, the adaptive grid and 
box-size are optimized |22|. When we vary the atomic 
mass from to/tohc — 1 to 0.692, the two-body binding 
energy for the KORONA potential changes by nearly 
ten orders of magnitude from E2 = —0.041 K to E2 = 
—8.6 x 10~ 12 K. The total ground state energies En for a 
given N near threshold, scaled by the atom mass m, are 
predicted to vary exponentially as a function of to |j| , 



ttiEn oc exp < 



Here, ml denotes the critical mass at which the N- 
body system becomes unbound, and af 1 ^ 1 parameters 
specific to the N particle system. The mass to is di- 
rectly proportional to the coupling strength K and in- 
versely proportional to the quantum parameter to which 
have been used previously to characterize LJ systems 
(see [I])). The functional form in Eq. (J2J) has been de- 
rived by Taylor-expanding the logarithmic derivative of 
the bound state wave function for small total ground 
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state energies about (to* 



(JV) 



1) [3 • Previous treatments 



neglected terms proportional to (m* — to) 1 with i > 1. 
To test whether this is justified, symbols in Fig. show 
the quantity — 1/ \n(m\E2\/mn c e) 1 where e denotes the 
well depth of the KORONA and LJ potential, respec- 
tively, as a function of m/mn e (e = 11.06 K for the KO- 
RONA potential and 10.22 K for the LJ potential). Di- 
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FIG. 1: Scaled ground state energies — 1 / \n(m\E2 \ /rriHet) for 
two 2D particles interacting through the KORONA potential 
(diamonds) and through a simple LJ potential (asterisks) as 
a function of m/muc- Solid lines show fits, which treat raf , 



and a 2 as parameters. Vertical solid lines indicate the 
mass ratio m* /mue at which the fit predicts the dimer to 
become unbound. The dotted line shows E2, Eq. (J^J, for 
the dimer interacting through the KORONA potential. The 
inset shows the ratios between the energies E\ and E2 (solid 
line), and between r e ff and a 3C (dashed line) for two particles 
interacting through the KORONA potential as a function of 



amonds show the ground state energies for the dimer in- 
teracting through the KORONA potential, and asterisks 
those for the dimer interacting through the LJ potential. 
Figure^shows clear deviations from a linear behavior, in- 
dicating that the term proportional to (m— ml ) 2 cannot 
be neglected for N = 2. Earlier studies [H El E3 did 
not see deviations from the linear behavior possibly be- 
cause the total ground state energies were i) varied over 
a smaller range, and ii) determined variationally. The 
qualitative behavior of the energies calculated using the 
KORONA and the LJ potential is similar, which implies 
that the non-linear behavior cannot be attributed to the 
difference in the long-range parametrization of the two 
two-body potentials [the LJ potential is for large r pro- 
portional to r~ 6 , while the KORONA potential contains 
additional terms proportional to r~ 3 , where j = 8 — 16 
(j even)]. 

Solid lines in Fig. ^ show fits of our scaled numeri- 
cal two-body ground state energies to Eq. J2J), treating 



,(2) J2) 



(2) 

and a 2 as fitting parameters and setting a 



(2) 



with i > 2 to zero. The fits predict that the two-body 
system interacting through the KORONA potential be- 
comes unbound at m* = 0.636(3)mnc and that inter- 
acting through the LJ potential at to* = 0.657(3)toh o - 
The numbers in brackets denote the uncertainties of the 
fit, which are obtained by including a varying number of 
data points in the fit. Tabled lists the fitting parameters 
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a\ ' and a 2 



(2) 



and their uncertainties. The critical 
masses m*~' for both interaction potentials are indicated 
in Fig. n by vertical solid lines. 

It is interesting to ask how well effective range theory 
describes the near-threshold behavior of the 2D dimcr. 
For a zero-range potential, the ground state energy E\ 
is determined by the 2D scattering length a sc (see, e.g., 
Ref. EfJ), 



E° 2 = - 



mat 



4exp(-2C), 



(3) 



where C denotes Euler's constant, C = 0.5772. Our def- 
inition of the scattering length a sc follows that adopted 
by Verhaar et al. [2^, that is, the scattering wave func- 
tion goes through zero at r = a sc . In the following we 
restrict ourselves to dimers interacting through the KO- 
RONA potential. The scattering length a sc , which we 
determine numerically, varies from 77.8 A for the 2D 
helium dimer to about 5.7 x 10 6 A for the most weakly- 
bound dimer considered with m — 0.692toh c - A dotted 
line in Fig. [2 shows E$, Eq. ©, using scaled dimension- 
less units. Near threshold, E% nearly coincides with the 
numerically determined energy E 2 (diamonds). At larger 
m, however, discrepancies are visible. To quantify these 
discrepancies, a solid line in the inset of Fig. ^ shows the 
energy ratio E2/E2, which varies from essentially 1 to 
0.8 for the atomic masses considered. Near threshold, 
the ground state energy E2 can be described to a very 
good approximation through a single atomic physics pa- 
rameter, i.e., the two-body scattering length a sc ; the im- 
portance of effective range corrections, however, increases 
with increasing mass. 

To account for a non-vanishing effective range, the 
right hand side of Eq. J3J has to be multiplied by 
exp(|S 2 | mr 2 s / (4h 2 )) |23,|2J|. The effective range r c g, 



h 2 



\E 2 \m 



In \Eo 



. exp(2C) 



h 2 



(4) 



can hence be evaluated if the scattering length a sc and 
the binding energy E 2 are known. For the dimers un- 
der study, r c g changes from 60 A for to/tohc = 1 to 
4900 A for to = 0.696toh c . For smaller masses the nu- 
merical accuracy of E 2 and a sc is not sufficient to reliably 
determine r e ff from Eq. Q . A dashed line in the inset of 
Fig. ^shows the ratio between r e g and a sc . Since the ef- 
fective range r e g increases with decreasing atomic mass, 
the energy scale associated with the effective range r e g, 
h 2 /(mr 2 s ) [25| . decreases with decreasing to (see below). 

We now turn to the study of 2D trimers interacting 
additively through the KORONA potential. Since the 
potential energy depends on relative coordinates only, we 
can separate off the center of mass motion. Restricting 
ourselves to states with vanishing total angular momen- 
tum reduces the number of degrees of freedom to three. 



The 2D Schrodinger equation can then be rewritten in 
terms of three hyperspherical coordinates. Here, we em- 
pl oy W hitten-Smith's democratic coordinates R, ■& and 
ip 26]. We determine the solution to the Schrodinger 
equation by first calculating a set of angular-dependent 
channel functions and then solving a set of coupled hy- 
perradial equations. Since our B-spline implementation 
closely follows that used in Ref. 27] for 3D trimers, we 
do not discuss it in detail here. We note, however, that 
the grandangular momentum operator in 2D differs from 
that in 3D (see, e.g., Ref. 28]) and that the range of 
the angular coordinate # changes from [0, tt/4] in 3D to 
[-tt/4, 7r/4] in 2D |2ll2iJ. We include up to twelve chan- 
nels in our calculations and estimate the uncertainty of 
the ground state energy E3, which, in addition to the 
number of channels included in the expansion, depends 
on the angular and radial grids employed and on the 
step size AR used to calculate the coupling matrix el- 
ements [22]], to be at the few percent level. 

Asterisks in Fig. Eta) show the ground state ener- 
gies I E3 1 on a log-scale as a function of the scaled mass 
to/toho- Since the system size increases with decreas- 
ing mass, the calculations become more involved as the 
mass decreases. For the masses considered, the total 
ground state energy varies over nearly four orders of 
magnitude from E3 = —0.180 K for to/tojjc = 1 to 
E 3 = -6.51 x 10~ 5 K for to/to Hc = 0.740. To test 
the applicability of Eq. asterisks in Fig. |2fb) show 
the scaled gound state energies — \n{m\E^\ / rriYi e e) as a 
function of the scaled mass to/tohc- As in the N = 2 
case, the scaled energies for N = 3 clearly show devi- 
ations from linear behavior. To determine the critical 
mass for the trimer system, we fit our data to Eq. (J2J 
treating to* 3 \ and ckp as fitting parameters. Ta- 
ble [fl summarizes the result of the fit, which is shown by 
a solid line in Fig. H£b). The uncertainties of the fit are, 
as in the dimer case, determined by including a varying 

number of data points in the fit. Notably, the predicted 

(3) 

critical mass to* ' for the trimer interacting additively 
through the KORONA potential nearly coincides with 
the critical mass to* for the dimer interacting through 
the KORONA potential. 

We now consider 2D systems with up to N — 7 parti- 
cles interacting additively through the KORONA poten- 
tial. For these larger systems, basis set expansion-type 
techniques become computationally unfeasible. We thus 
solve the 2D many-body Schrodinger equation by alter- 
native means using the essentially exact DMC technique 
with importance sampling [30j. The DMC technique al- 
lows ground state energies and structural properties to 
be determined. Since the numerical solution of the time- 
independent Schrodinger equation is based on a stochas- 
tic process, expectation values can only be determined 
within a statistical uncertainty. This statistical uncer- 
tainty can be reduced by increasing the computational 
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FIG. 2: (a) Symbols show the total ground state energies En 
for N — 2 (pluses), TV = 3 (asterisks), TV = 4 (diamonds), 
TV = 5 (triangles), AT = 6 (squares), and N — 7 (crosses) 
on a logarithmic scale as a function of m/mHe for the KO- 
RONA potential. Dotted lines connect data points for fixed 
TV to guide the eye. The upper solid and dashed lines show 



h /(mr e ) for r e 



>'cH' 



and 



r v dw, respectively. 



lower solid and dashed lines show O.lh /(mr e ) for r e 



The 

= r cB 

and r e = r v( jw, respectively. A vertical solid line indicates 
the mass ratio m* /ffiEe = 0.636, at which our fits predict 
the 2D dimer to become unbound; the critical masses for the 
larger systems nearly coincide with that for the dimer. (b) 
Scaled total ground state energies — 1/ ln(ro|-Ejv|/mHee) as a 
function of m/muc using the same symbols as in (a). Solid 



lines show fits (treating mi , and a y - y ' as parameters 

for TV = 2 — 5, and treating m* and as parameters for 
TV = 6 - 7). 
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efforts. Our DMC energies for the 2D helium trimer 
and tetramer agree with those calculated by Vranjes and 
Kilic |3l| . Furthermore, our DMC energies for the trimer 
agree to within the statistical uncertainty with those cal- 
culated by the hyperspherical B-spline treatment (see 
above). In the following, we report the ground state en- 
ergies En for TV > 3 calculated by the DMC method as a 
function of the atomic mass. As we approach the thresh- 
old regime, the DMC calculations become more difficult 
since the kinetic and potential energy nearly cancel. 

Symbols in Fig. show the DMC ground state en- 
ergies En as a function of the mass ratio m/mHe for 
TV = 4 — 7. Panel (a) shows the energies on a loga- 
rithmic scale (to guide the eye dotted lines connect data 
points for the same TV). To investigate how well the 
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FIG. 3: Energy ratios En+i/En as a function of m/miie 
for TV = 2 (uppermost curve) through TV = 6 (lowermost 
curve) calculated using the KORONA potential. Errorbars re- 
flect the statistical uncertainties of our DMC energies En for 
TV > 4. Dashed horizontal lines on the left hand side indicate 
the energy ratios for zero-range interactions: E3/E2 = 16.52, 
Ei/Ei = 11.94 and E S N+1 /E% = 8.57. A vertical solid line 
indicates the mass ratio mi 2 ^ /mnc = 0.636, at which our fits 
predict the 2D dimer to become unbound. The critical masses 
for larger systems nearly coincide with that for the dimer. 



functional form given in Eq. @ applies to systems with 
TV > 3, panel (b) shows the scaled DMC ground state en- 
ergies — 1/ ln(m|75/v|/mHee)- Statistical uncertainties of 
the DMC energies (not shown) are smaller than the sym- 
bol size. Deviations from a linear behavior are, although 
less pronounced than for the dimer and trimer, visible for 
TV = 4 and 5; consequently, we fit our DMC energies to 



Eq. J2J> treating 



(TV) (AT) , (AT) „ . 

i* , a\ and cir, as fitting parame- 
ters. For the clusters with TV = 6 and 7, our scaled DMC 
energies depend to a good approximation linearly on the 
mass; hence we use only two fitting parameters, mi N ^ 



and 



AN) 



The fitting parameters are summarized in Ta- 



ble |U and the fits are shown by solid lines in Fig. HJTb). 
The statistical uncertainty of the fits is, as in the dimer 
and trimer case, estimated by including a varying number 
of data points in the fit. We speculate that, if we were 
able to obtain accurate DMC energies closer to threshold, 
non-linear behavior of the scaled ground state energies for 
TV > 6 would be, similarly as for the dimer and trimer, 
visible. The fits for TV = 4 — 7 predict critical masses 
toI^/tohg between 0.62 and 0.63, which agree within 
their uncertainties with those predicted for TV = 2 and 3. 
Our analysis thus confirms that the critical mass mi at 
which 2D systems, interacting additively through a given 
two-body potential, become unbound is the same for all 
system sizes [l(J, |llj, |l2| . 

To further investigate the near-threshold behavior for 
TV > 2, symbols in Fig. |3| show the ratio between the 
ground state energies for systems with TV+l and TV atoms 
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as a function of m/mn e - To guide the eye, dotted lines 
connect data points for Em+i/E^ with the same N + 1 
and N but different to. Errorbars, which increase with 
decreasing mass, reflect the statistical uncertainty of our 
DMC energies for N > 4. At large to, the energy ra- 
tios approach a constant, in agreement with earlier stud- 
ies |l2|. Since I-E2I decreases with decreasing mass while 
the van der Waals length r v( jw changes only little [see 
dashed lines in Fig. Eta)] , dimers with small m effectively 
have a shorter range than those with large m. Accord- 
ingly, systems with small m should be better described by 
zero-range models than those with large to. Horizontal 
dashed lines on the left hand side of Fig.EJindicate the en- 
ergy ratios El/Ei = 16.52 0, E s jEl = 11.94 [l| and 
E s N+l /E s N = 8.57 H3 for 2D systems interacting through 
additive zero-range potentials. While the energy ratio 
E3/E2 calculated for the helium- like few-body systems 
with small to is close to the value predicted by the zero- 
range treatment, none of the energy ratios En+i/En for 
N > 3 is. Although the energy ratios for N > 3 in- 
crease with decreasing mass, it is not clear whether they 
approach the values predicted for the zero-range model. 

Recall that zero-range treatments should be applica- 
ble when the energy scale associated with r e is much 
larger than \En\, where r e is either given by the van 
der Waals length r v aw or by the effective range r e g. 
Solid and dashed lines in panel (a) of Fig. El show 
the energy h 2 /(mr 2 ) for r e = r c g (upper solid line) 
and for r e = r v( jw (upper dashed line), respectively. 
Since \En\ <C h 2 /(mrl) for zero-range treatments to be 
valid, the lower solid and dashed lines show the quan- 
tity 0.lh 2 /(mr 2 ) for r e — r e g and r e = r v( jw, respec- 
tively. Figure Eta) indicates that, of the systems consid- 
ered, only those with N = 2 and 3 have total ground 
state energies \En\ that are smaller than the lower solid 
and dashed lines. Although the dimer near threshold is 
well described by a zero-range treatment (see Fig.^), the 
corresponding systems with N > 3 cannot be properly 
described through simple contact interactions. This is 
due to the fact that the total ground state energies vary 
exponentially as a function of N and that the effective 
range r c g increases with decreasing to. To check whether 
the latter is specific to the KORONA potential studied 
here, we also calculated the effective range r e ff for the 
LJ potential discussed above and for a realistic tritium- 
tritium b 3 S+ potential. For the three potentials consid- 
ered, the effective range behaves similarly as a function of 
the two-body binding energy. Our findings suggest that 
the regime where the A-body total ground state energies 
Ejf, A ;§> 2, can be properly described by zero-range 
models might, in general, be hard to reach for 2D van 
der Waals systems. 

Finally we remark on an earlier variational study Q, 
which concluded that artificial "bosonic helium 3" clus- 
ters in 2D become bound for a minimum of about twelve 
atoms. Although the figures in this paper contain ener- 
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mi N) a[ N) 


(N) 


LJ 


2 


0.657(3) 0.66(2) 


-0.53(6) 


KORONA 


2 


0.636(3) 0.68(2) 


-0.55(6) 




3 


0.637(8) 0.86(5) 


-0.53(14) 




4 


0.62(2) 0.89(8) 


-0.25(12) 




5 


0.62(2) 1.06(10) 


-0.26(20) 




6 


0.62(3) 1.14(10) 






7 


0.63(2) 1.30(7) 





TABLE I: Fitting parameters for 2D clusters interacting 
through the Lenard- Jones potential (for N = 2 only) and 
through the KORONA potential (N = 2 — 7), respectively. 
We treat three fitting parameters for N = 2 — 5, and two 
fitting parameters for TV = 6 — 7. The numbers in brackets 
give the uncertainties of the fitting parameters, which are ob- 
tained by including a varying number of data points in the 
fit. 



gies for to < m( 3 He) for only a few system sizes, our cal- 
culations show that bosonic 3 He clusters are self-bound 
for all A considered, in agreement with Ref. 31]. (We 
excluded some energies for small m from the figures and 
the fits since their statistical uncertainties are very large.) 
Furthermore, the existence of a universal critical mass 
indicates that "bosonic helium 3" 2D clusters, which are 
often studied to complement the treatment of fcrmionic 
helium 3 systems, are bound for all A. Our calcula- 
tions emphasize that great care has to be taken when 
variational calculations are used to investigate the near- 
threshold regime (see also Ref. 0). 
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